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Abstract — As computing power and algorithmic 
advances have evolved significantly in the recent past, it is 
now feasible to solve complex electromagnetic (EM) 
problems involving scattering, radar cross section, 
antenna design, microwave circuit design, artificial EM 
materials, etc. using full-wave numerical methods. Several 
general-purpose commercial software packages are 
routinely used in industry in all these domains for EM 
analysis or design. However, the task of processing large 
sets of data output from these design studies and analyses 
is generally beyond the realm of commercial software 
packages, and the designer spends many hours writing 
problem-specific computer programs to extract the 
desired performance parameters. Determination of 
coupling coefficients or the unloaded quality factor of a 
dielectric resonator, de-embedding feed lines from 
antenna currents, removal of discontinuity effects, the 
extraction of equivalent circuit models, etc. are some 
examples where auxiliary processing is needed for the 
extraction of EM parameters of interest. The same 
considerations apply to the parametric analysis of 
measured data in the presence of noise. This paper 
presents a versatile data-driven spectral model derived 
from a state-space system representation of the computed 
or measured EM fields, from which all the parameters of 
interest can be extracted. Several examples will be 
presented to demonstrate the usefulness of the proposed 
approach for parametric extraction in EM problems. 


Index Terms —Electromagnetic Scattering, Finite Difference 
Time Domain, Finite Integration Time Domain, Microwave 
Circuit, Microwave Resonator, Signal Extrapolation, Spectral 
Estimation, State Space Method, System Identification, 
Transient Analysis, Transmission Line 


I. INTRODUCTION 


With prolific algorithmic advances and significant 
improvement of computational resources in recent years, full- 
wave electromagnetic (EM) simulation techniques, such as the 
method of moments (MoM) [1], the finite element method 
(FEM) [2], the finite-difference time-domain (FDTD) method 


[3] and the transmission line matrix method [4], are 
increasingly used to solve complex EM problems involving 
scattering, radar cross section (RCS), antenna design, 
microwave circuit design, artificial EM materials, etc. These 
rigorous techniques account for physical phenomena such as 
surface-wave coupling and radiation, dispersion, frequency- 
dependent metallization and dielectric losses, proximity 
effects and near-field coupling. Although EM simulation 
methods are computationally intensive in practice, large-scale 
GPU and distributed computing enables commercially 
available full-wave EM analysis software packages to be 
applied to challenging problems such as scattering by 
electrically large objects (e.g., computing the RCS of an 
aircraft) and modeling the EMI between interconnects on a 
printed circuit board containing dense RF and/or high-speed 
digital circuits. The outputs of these software packages are 
generally voltages, currents, scattering (S) parameters or some 
derived EM field variable such as the RCS. In design studies 
large datasets of these output parameters need to be processed, 
for example, to remove the effects of transmission line 
discontinuities on the S-parameters or to extract the modal 
response of a particular wave species and isolate the scattering 
centers causing undesirable RCS. Therefore, it is very 
desirable to investigate independent signal models (in time 
domain) and spectral models (in frequency domain), which 
accurately represent the simulated or measured output fields. 
One can then utilize this model for removing the effects of 
discontinuities or parasitic scattering mechanisms so that the 
desired EM performance can be characterized accurately. 


Several researchers have extracted parameters of interest 
from scattered fields or circuit response using signal 
processing techniques such as Prony’s method [5], [6], pencil 
of functions [7]-[11], autoregressive moving average 
(ARMA) [12], [13], estimation of signal parameters via 
rotational invariance techniques (ESPRIT) [14]-[16], multiple 
signal classification (MUSIC) [17], and the state space 
method (SSM) [18]-[21]. Applications include computation of 
complex natural resonances and eigenmodes [22]-[28], 
impulse response characterization of time-domain signatures 
[29]-[34], broadband circuit parameter extraction [35], [36], 
identification of radar target’s features [37]-[40], extraction of 
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biomedical vital signs from UWB radar measurements [41 ]- 
[43], and location of buried targets using ground penetrating 
radar [44]. The basis behind such signal processing 
application is that the EM field scattered by an object can be 
adequately represented as a sum of damped sinusoids, whose 
amplitude and phase are closely related to the physical 
parameters of interest. For example, poles of the transfer 
function for the exponential signal model locate discrete 
scattering centers useful in object typing and feature 
identification in radar target identification. 


This paper presents a versatile data-driven spectral model 
derived from a state-space system representation of the 
computed or measured scattering parameters and EM fields, 
from which all the parameters of interest can be extracted. 
Parametric extraction from measured data is especially 
challenging because of noise and random measurement errors. 
The efficacy of spectral estimation methods to extract 
parameters of damped sinusoids embedded in noise has been 
studied by many researchers. For example, performance 
analysis of MUSIC is treated in [45], [46], estimation of the 
direction of arrival of radar returns using ESPRIT is presented 
in [47], a performance study of matrix pencil method in the 
presence of noise is described in [48], [49], and sensitivity 
analysis of the state space method is treated in [50], [51]. SSM 
[18]-[21] has been extended to ultra-wideband coherent 
processing of range-Doppler data for radar target 
identification and validated with static range measurements 
[52]-[54]. A comparative analysis between state space and 
matrix pencil methods shows similar performance with 
comparable accuracy when implemented on harmonic 
retrieval problems with noisy data [49], [55], [56]. 


The reader is referred to [40]-[43], [52]-[54] for the 
application of SSM to feature identification using monostatic 
RCS measured data on stationary targets as well as human 
subjects. As we have shown in [40] using measured data on a 
large conical target with and without a lossy dielectric coating, 
the isolation of electromagnetic wave species of interest, such 
as creeping waves, multiply diffracted waves and specular 
scattering, yields a better understanding of the physics behind 
wave propagation around curved dielectric or coated 
structures, thereby improving the accuracy of feature 
extraction or target identification. For example, it has been 
shown that the coating enhances the creeping wave 
contribution significantly in certain directions compared to the 
metallic cone, suggesting its efficacy in absorbing EM waves 
[40]. The biomedical UWB radar [41]-[43] uses narrow 
pulses to probe the human body and detect tiny 
cardiopulmonary movements by spectral analysis of the 
backscattered EM field. With the help of super-resolution 
spectral algorithms [52]-[54], UWB radar is capable of 
increased accuracy for estimating vital signs such as heart and 
respiration rates in adverse signal-to-noise conditions. A 
major challenge for biomedical radar systems is detecting the 
heartbeat of a subject with high accuracy, because of minute 
thorax motion (less than 0.5 mm) caused by heartbeat. The 
problem becomes compounded by EM clutter and noise in the 
environment. We have shown that SSM processing of the 
UWB radar data on a stationary human subject consistently 
produces accurate estimates of the vital signs for several 


channels of UWB data without producing harmonics and 
inter-modulation products that plague signal resolution in 
widely used FFT spectrograms [41]. 


The state space method can be applied in either time domain 
or frequency domain. In this paper, we focus on the frequency 
domain problems. The reader is referred to [34] for the 
application of SSM to time-domain EM problems. An 
attractive feature of SSM 1s its ability to identify and associate 
a small number of poles of the system transfer function with a 
specific scattering mechanism or modal response, such as 
scattering parameters at the transition between a coaxial 
connector and a microstrip transmission line in measured data, 
and isolated scattering by edges and seams from the composite 
RCS of a large body. Thus, the desired field parameter can be 
extracted or estimated from synthetic or measured data using 
a linear system with relatively small model order. Illustrative 
examples will be presented to demonstrate the usefulness of 
the state space method for parametric extraction in EM 
problems. We describe the problem formulation and the 
modeling process in detail so that an uninitiated reader can 
follow the steps to program and execute the algorithm and 
generate results of interest. 


In the first example, we consider a planar dielectric slab 
illuminated by a plane wave at normal incidence and isolate 
the reflection off the front face using range-classified poles 
pertinent to the specular reflection. This process is akin to de- 
embedding transmission lines at the ports to evaluate the 
circuit behavior of an antenna or a discontinuity. In the second 
example, appealing to the canonical problem of Mie scattering 
by a sphere, creeping waves are extracted using SSM to model 
the RCS response, and validated against high-frequency 
asymptotic approximations. Next, random white Gaussian 
noise is added to the Mie series solution, and Monte Carlo 
simulation is performed to examine robustness of the SSM 
estimates to noise. Numerical considerations such as dynamic 
range, signal-to-noise ratio (SNR) and model order 
determination are addressed in detail for both examples. 


The paper is organized as follows. In Section II the state 
space algorithm is reviewed in detail following [52], [53]. 
Section III presents basic numerical considerations on the 
state space approach, such as estimating the SNR and the 
model order. Section IV presents illustrative examples on 
parametric extraction for EM problems, as discussed earlier. 
Factors affecting accuracy and numerical efficiency of the 
proposed technique are discussed. Finally, some concluding 
remarks are provided in Section V. 


Il. STATE SPACE METHOD 


In recent years, there has been a great deal of attention 
devoted to model-based eigen-decomposition methods 
derived from a state-space realization of the system 
identification problem, which is defined as the determination 
of the internal states of a linear time-invariant (LTI) system 
given a Set of inputs and outputs [18], [21], [52]-[62]. The idea 
of using state-space methods to estimate frequencies and 
amplitudes of damped sinusoids was first suggested in [18], 
where Kung et al. developed a system identification approach 
based on singular value decomposition (SVD) for the 
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harmonic retrieval (or spectral estimation) problem. The 
foundation of state-space signal modeling is based on 
representation of a linear rational system, popular in linear 
systems and control theory [63]-[65], 1n which the difference 
equations for the discrete-time (or discrete-frequency) signal 
are converted into state equations and the model parameters 
are estimated in terms of the state matrices characterizing the 
system. State-space parameterization enables reduction of 
parameter sensitivity, as demonstrated by several examples 
relevant to the sinusoid retrieval problem in the tutorial article 
by Rao and Arun [21]. Unlike polynomial-based signal 
processing methods such as Prony’s, MUSIC and maximum 
likelihood estimation, in which the frequencies are obtained 
by cumbersome root-search of the polynomials, state-space 
methods simultaneously yield complex amplitudes (with 
initial phases) and frequencies directly from three state 
matrices. The frequencies are calculated from the eigenvalues 
of the state transition matrix, and the amplitudes are derived 
by modal eigen-decomposition of the state equations using 
two auxiliary matrices, the control and observation matrices. 





As shown in the sequel, because the decomposition of the 
data into a sum of damped sinusoids gives rise to an LTI 
system, the state matrices can be readily computed via a low- 
rank truncation of the Hankel (or forward-prediction) matrix 
representing the data. Thus, each entry of the Hankel matrix 
can be expressed in terms of the impulse response derived 
from the state space matrices. The SVD of the Hankel matrix 
is a product of the observability and controllability matrices 
which yield all the model parameters in the harmonic retrieval 
using state dimensions corresponding to the number of 
scatterers embedded in the data. 


A. ARMA Signal Model 


The scattered field output data sequence y(k) comprises N 
uniformly spaced frequency samples (see (1) below), each 
represented as a sum of M complex sinusoids (or scattering 
centers) corrupted by measurement noise w(k), assumed to be 
white Gaussian with zero mean. In deterministic data 
modeling we assume that w(k)=0 and characterize the 


impulse response, thus estimating the system parameters using 
only the output. When the input is known and needs to be 
considered (e.g., short pulses or modulated waveforms) the 
system can be identified using both input w(k) and output y(‘) 
[34]. For completeness we retain w(h) in the formulation even 
if we may not consider it in a given case. Thus, over a given 
bandwidth, the signal measurements at N frequencies are 
modeled as 


M 
PER ae PEN ews kale Ny 
i=l 
Y(k) = yk) + w(k). 
The difference between the “true” signal (or truth), y(‘), and 
V(k), 


measurement or modeling error, w(k). The model parameters 


the state space model, is the random noise, 
a, and a, represent the amplitude and its spectral rate of 
decay (for negative a,) or growth (for positive a), 
respectively, for the i-th scattering center. Both positive and 


negative a, are needed to model the peaks and dips in the 


frequency response (for example at resonance). It is 
emphasized that positive œ, does not imply instability or 


violation of any physical constraint such as the radiation 
condition. The parameter r, denotes time delay of the i-th 


scatterer at the observation point and is related to range R, by 
T, =2R./c, where c is the speed of light. Equation (1) is set 


up in terms of the monostatic RCS that an observer would 
measure for a given target at range R. Therefore, 2R, is 


simply the round-trip distance between the transmitter and the 
target. If bistatic radar is used, then 2R, should be replaced 


by R. +R, 
from the transmitter and the receiver, respectively. It follows 
that 277, f, =(@,/c)\(2R,) = 8,(2R,) represents the phase of 


the i-th sinusoid, with Ø, being the phase constant at the 


where R, and R, denote the distance to the target 


frequency f,. The frequency vector is specified in terms of 


the carrier frequency f, as 


fe=f +&-DAS, k=1,--,N, (2) 
where Af is the sampling frequency. The primary interest in 


state space system identification 1s to estimate the parameters 
a,,a@,andR, (or z,), which are embedded in the data 
sequence j(k). The latter two parameters are computed from 
the eigenvalues of an open-loop state matrix, to be defined 
shortly. Once these parameters are estimated, the amplitudes 
a, can be readily derived from the state matrices using a 


modal decomposition method based on least squares [53]. 
B. ARMA Transfer Function 


The ARMA model in (1) is interpreted as an LTI system 
with input given by the sequence w(k) and the output by the 
data sequence y(k). The goal in system identification is to 
compute the coefficients of the ARMA transfer function (TF) 
characterizing the discrete-frequency signal model in (1). 
Taking the z-transform of (1) after substituting f, from (2), 


we obtain 
a B. 
Y(z)= >) ————— +W (z 3 
@= m a TO (3) 
A 4, 
B, Sge pih a Ag” (4) 
P; =Q, + J2AT, (5) 


In (3) - (5) B, denote the complex amplitudes, 0, is the initial 
(static) phase, and p, are the poles in the complex z-plane. 


The summation in (3) results in a TF comprising rational 
polynomials in the numerator and denominator: 


M-1 
c+) CZ 
yg de 


a = M 
W(2) l- >. d z' 
i=] 


N(z) 
D(z) 


I> 


T(z) 





(6) 
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The TF 7(z)has M poles and M —1 zeros, located at the 
roots of D(z) and N(z), respectively. The TF in (6) represents 


a special case of a more general ARMA TF given by 





Q 
+ oz" 
či Lo RA 
e pe (7) 
1- > azn á 
i=] 
which has P poles and Q zeros (P and QO can assume arbitrary 
integer values). The special case Q = P—1 listed in (6) is also 
known as “strictly proper” TF and the Prony’s model belongs 
to this category. Another important special case is the purely 
Auto Regressive (AR) TF with Q =0. 
The input-output relationship for the general ARMA model 
with the TF in (7) is characterized by the difference equation 


P Q 
y(k) = 2,4 y(k —i)+ Ddcjm(k -j)+cow(k) (8) 


Eq. (8) is valid for the TF in (6) as well with P =M and 
Q = (M —1). From linear systems and control theory [63], 
[64] one can show that the difference equation in (8) may be 
written alternatively in terms of the state-space description 
characterized by the difference equations 

x(k +1) = Ax(k)+ Bw(k) (9) 
Vk) = Cx(k) + Wk), (10) 
is the state, AEC” is the state 


transition matrix, BeC””' and CeC”” are constant 
matrices known as control matrix and observation matrix, 
respectively [64]. Our goal is to identify the matrices A, B and 
C given the data sequence (4) and the input w(k). The transfer 
function, 7(z), is obtained by taking the z-transform of (9) and 
(10) and evaluating the ratio Y(z)/W(z): 


where x(k) e C*“ 


T(z)=C(zI-Ą `B +1. (11) 


It follows from (11) that the poles of the model (1.e., the roots 
of D(z) in (7)) are the eigenvalues of the open-loop matrix A 
and the zeros (1.e., the roots of N(z) in (7)) are the eigenvalues 
of the matrix (A— BC) [52], [53]. 


C. State Space System Identification 


The transfer function in (7) can be written equivalently in 
terms of its impulse response sequence as 


T(z) = A(O) +A +-+ A(n)” +- (12) 


which is aptly referred to as the infinite impulse response (IIR) 
transfer function. The expansion of the inverse matrix in (11) 
into an infinite series yields 


(zI -A` = Iz’ + Az”? +A’ z” + (13) 
By inserting (13) into (11) we obtain 
T(z) =1+CBz'+CABz* +CA°Bz? +» (4 


Equating coefficients of like-powers of z in (12) and (14) and 
realizing that h(k) = y(k) for the impulse input, we obtain 


y(0)=1 

yQ) = CB 

y(2) f CAB (15) 
y(k) =CA*'B 


Therefore, the relationship between the impulse response of 
the model and the state-space parameters for any positive 
value of k is defined by 


y(k) = CAB, k >Q. (16) 


Eq. (16) indicates that a Hankel (or forward-prediction) 
matrix, H, formed from the IIR sequence of a system as 
ya) y2) y6) 
2 3 4) «+ 
pee) POr (17) 
yB) IH I5) 


can be decomposed into a product of two matrices given by 


C 
CA : F 
Maen [B AB &B -—|£0rF, (18) 
where Q and I are known as observability and 


controllability matrices, respectively [65]. It 1s important to 
note that despite the infinite dimensions of H in (17), in 
practice the impulse response is always finite. Thus, for a 
given set of measurements the rank of the Hankel matrix H, 
and by inference the rank of Q and T, will always be finite. 


As described later in this Section, Q and IT , consequently H, 
can be truncated to low-rank matrices with rank r< M, 


where M is the number of complex sinusoids in the model (see 
(1)). In summary, the Hankel matrix H may be interpreted as 
an operator constructed from a set of measurements y(k) that 
maps the past input vector w to the future output y*. Causality 
of SSM is explicit in (8) - (10), which emphasize that the 
current output is dependent on the past output and the current 
as well as the past input signals. Next, we present a method 
to derive the state-space matrices from the Hankel matrix 
constructed using a finite set of output data samples. 


The first step in computing the triplet (4, B, C) is to form 
the Hankel matrix using the available data samples 
y(k), A =1,2,---N. 


yQ) y(2) y(L) 
ee ra n MeT 1) (19) 
y(N-L+1) y(N-L+2) y(N) 


where the parameter L denotes length of the correlation 
window, heuristically chosen to be L=[N/2], and the 
brackets denote the smallest integer less than or equal to the 
inserted quantity. Note that the Hankel matrix in (19) is a 
truncated version of the IIR in (17). Subspace decomposition 
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methods exploit the eigenstructure of Hankel matrices to 
compute the state matrices of the LTI system and estimate the 
signal models [18], [21], [25], [26], [52]-[56], [60]-[62]. 
Accordingly, we partition the Hankel matrix into signal and 


noise subspaces using the singular value decomposition of H. 
Organizing the SVD in terms of the singular values fo} in 


l 


descending order of magnitude, the Hankel matrix may be 


written as 
> 0 NV" è 
H =|U, A [7 -v27 (20) 
0 LTV, 


where the subscripts s and n denote the signal and noise 
subspaces, respectively, and the asterisk refers to matrix 
conjugate transpose. The matrices U, and U, are the signal 


and noise components, respectively, of the left-unitary matrix 
[U, U, |. Likewise, V and V, denote the signal and noise 
of the 
|V, y|. Furthermore, 2, and X, are diagonal matrices 


components, respectively, right-unitary matrix 


comprising the signal and noise singular values, respectively. 
It is understood that the signal components in the SVD are 
entirely characterized by the dominant singular values È. 
The classification between signal and noise subspaces is 
achieved by parsing the singular value spectrum in descending 
order of magnitude and removing the noise singular values ÈX, 
[18], [50]-[53]. To increase the accuracy of the state-space 
matrices, the Hankel matrix H may then be truncated by 
suppressing the noise singular values and their associated 
unitary matrix components. This results in a reduced-rank 
approximation to the Hankel matrix in (20), obtained by 
retaining only the dominant singular values (cf. [51]), 1.e., 


AaU zy. (21) 


Let the computed singular values in X be arranged in 
descending order of magnitude 


(22) 
where r is chosen subject to the threshold [53] o,/o, ~ 10°”, 


O,>0,> + >a, >0 


and p is the number of significant decimal digits in the truth 


data. For example, if the data is accurate to three significant 
digits, then the singular values with upper bound r for which 
the ratio is less than 0.001 are considered as noise singular 
values and excluded from the model. The largest index r of the 
singular values in the signal subspace is the rank of the 
Hankel matrix. It can be shown that the largest retained 
singular value minimizes the error between the Hankel 
matrices H and H in the spectral norm sense [49], i.e., 


O, ~|H -Ë 








(23) 


where the subscript s denotes the spectral or Lz norm [64]. As 
demonstrated in Section II, magnitude of the dominant 
singular values in the signal subspace may be conveniently 
used to estimate the “optimal” model order M in (1). Next, we 
address how one can compute the state matrices from the 
SVD. 


Akin to (18) H is obtained in factored form as 


9 
S 


HeU Z =r (24) 
By using the balanced coordinate transformation method 
proposed in [18], one can compute the finite-rank 


observability matrix Q and controllability matrix I from 
the SVD in (24). These matrices are given by 


Õ=U IŻ and Õ =X% 


Cus 


(25) 


Then, the open-loop matrix Ae can be derived from 


either Q or I as shown in the Appendix. If the derivation of 
A is based on the observability matrix Q , then [40] 


~ ~ = wes A 
A= (aO) aa (26) 


The matrices Q_, and Q_,, in (26) are obtained by deleting 


the first and last rows, respectively, of Q. Alternatively, in 
terms of the controllability matrix, A is given by [40] 


A= als (I al 


= —cl 


(27) 


~ 


where the matrices [_., and I'_,, are obtained by deleting the 


-rl 
first and last rows, respectively, of I. The reader is referred 
to the Appendix for the least-squares computation of B and C, 
also using either Q or I. As noted earlier, the matrices B and 
C contribute to zeros of the transfer function. In the AR SSM 
that we have employed in [39], these two matrices are not 
used, because the system is entirely characterized by the poles 
determined from A. 


D. Modal Decomposition 


Once the state matrices (4,B,C) are known, the model 
parameters in (1) may be computed using a modal (or eigen-) 
decomposition method [52], [53]. If the complex eigenvalues 
of A are assumed to be distinct, one has 


At A={A, A, At (28) 
The magnitude of the eigenvalues {A} determines amplitude 
decay/growth rate (with respect to frequency) of the system 
output, and their phase provides the time delay r, or the range 
R, in (1). As noted in (11), the eigenvalues of A represent the 
poles of the transfer function, T(z). The eigenvalues {A} and 


the eigenvectors ly, | of the state matrix A satisfy 

Ay, =Ay,; i=1--,M (29) 
Now let us form a modal matrix Y, from these eigenvectors 
by stacking them row-wise: 


Y = ly, Wa Wir | (30) 
It follows from (28)-(30) that 
Aa: (31) 
A=diag{ A, Aaj l, (32) 
Therefore, we may calculate A from (31) as 
A=P AP (33) 


The entries on its main diagonal are exactly the eigenvalues 
of the state transition matrix A. Therefore, the modal 
decomposition approach provides valid information to 
characterize the target, and it can be used to identify point 
scatterers embedded in the data set. The IIR approximation, 
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y(k), of the data sequence y(k) can be computed in terms of 
the state matrices using (16) as 


$(k) = CAB, kb =1,2,--N. (34) 


After substituting for A from (31) and realizing that 
A” = PA", we obtain 
P(k) = CV AL PB. (35) 


The amplitudes of the scatterers in (1) may be computed from 
(34). Let 


Psj v vy] (36) 
where the superscript T denotes matrix transpose. After 
inserting (30), (32) and (36) into (34), we obtain 

M 
Pk) =X (Cy, OBA T, k=l N (37) 
i=l 


Using (1), we deduce that the magnitude and phase of the 
eigenvalues (poles) are related to the model parameters a, 
T,, and R., respectively, by 


log|A, | 
æ, = el E , a pa , 
Af 2mAf AnAf 





=1,---,M 


(38) 
In (38), ø refers to phase of the eigenvalue 4. The 
amplitudes a, defined in (1) are then obtained from (37) using 
these eigenvalues, the corresponding eigenvectors, and the 
state matrices B and C: 
Cy,)(v.B 
gal aes F 
CaM 
The frequency dependence of the amplitudes and their 
decay/growth rate is emphasized in (38) and (39). 


iste M. (39) 


Next, we discuss the physical interpretation of the state 
Space parameters. The correspondence between the scattered 
response and its state space representation in terms of complex 
sinusoids becomes apparent when one compares (1) and (37). 
The poles A, of the transfer function yield the localized 


scattering centers on the scattering object. As derived in the 
Appendix, the matrices B and C follow from least squares 
solution applied to linear systems represented by the 
observability and controllability matrix, respectively (see (63) 
and (67)). Therefore, the matrices B and C are directly related 
to the observability and controllability of the linear system in 


(34). The observability matrix Q affects the shape of the 
natural response of the scattered field (directly proportional to 
the natural modes or eigenvectors of the state transition matrix 
A) within a given frequency band. The controllability matrix 


I limits the effect of an input such as an impulse on the 
scattering response. In other words, the controllability matrix 
determines how much of the excitation is coupled to the 
natural modes of a particular scattering center within the 
defined bandwidth. 


III. NUMERICAL CONSIDERATIONS 
A. Estimation of the Model Order 


In any spectral estimation problem, selection of model order 
of the underlying system is a critical issue. The difference 
between the desired and modeled signal, as defined in (1), is a 
combination of model mismatch error, measurement error, 
and noise. If this error is Gaussian distributed, then the 
minimization of the error criterion in (23) represents the 
maximum likelihood estimate (MLE) for the ARMA problem 
[18], [66]. In case of MLE, the Akaike information criterion 
(AIC) or minimum description length (MDL) criterion are 
often used for model order determination [67]-[69]. An 
estimator that yields the true number of signals with 
probability one as the sample size increases to infinity is said 
to be “consistent.” It has been shown in [68] that MDL is a 
consistent estimator of order, whereas AIC 1s inconsistent and 
often overestimates the model order. 


SSM has been used with simulated as well as measured 
data. In deterministic data modeling such as simulated data 
where the SNR is quite high, we employ the singular value 
matrix xX to estimate model order in the SSM. The large 
singular values in correspond to strong signal components, 
while the smaller values are generally attributed to noise. For 
low noise levels, there is a sharp transition between the large 
and small singular values. This transition point can be used as 
an estimate of the model order. At higher noise levels, 
especially for measured data, the transition from large to small 
singular values may not be well-defined, making model order 
estimation more difficult. In this case, probabilistic methods 
such as AIC and MDL criteria have been used for model order 
estimation (cf. [33]). No general criterion exists to determine 
which of these methods gives the “optimal” model order, and 
none of them provide necessary and sufficient conditions to 
guarantee that all the required signals are included in the 
model. It is expedient to consider the underlying physics of 
the scattering mechanism to guide the determination of the 
number of scatterers using a combination of these methods, 
and to examine the goodness of fit for each model order by 
computing the mean square error in the model [40]. For 
convenience, the equations to compute AIC and MDL from 
the singular values, o,, of the Hankel matrix approximation, 


H, are summarized below [68]. 


M iis 
—(M-r 


~ 





AIC(r) = —2(M —r)N In a 7 + 27r(2M —-r) 
Do 
MT ey 
(40) 


IT Gi 
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i=r+1 
1 M 


MDL(r) =—-(M -r)N In + - r(2M -r)InN 





O; 
Mak Sa 


(41) 
In the above, M is the maximum model order in (1) and N 


is the number of data samples. For a given number r of 
signals, the term in the parentheses simply denotes the ratio of 
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geometric mean to arithmetic mean of the smallest M -r 
singular values, and it is a measure of the error between the 
truth and the model. The number of signals, 1.e., the model 
order, is determined as the value of re {I,2,---(W—1)} for 


which either the AIC or the MDL is minimized [68]. 


iü Singular Value Spectrum vs. Model Order 
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Fig. 1. Singular value spectrum vs. model order for measured RCS 
on a cone. Circles denote MDL and diamonds denote AIC [40]. 


As an example, Fig. 1 shows the singular value spectrum 
(amplitudes in descending order of magnitude) plotted against 
the model order for experimental RCS data on a right circular 
conical stationary metallic target with and without a lossy 
dielectric coating [40]. It is seen that transition from strong 
signal components to small singular values is not well-defined 
because of measurement noise. AIC yields an order of 15 and 
13 for bare (metallic) and dielectric-coated cones, while MDL 
yields model orders of 4 and 7, respectively. Due to 
measurement noise, AIC usually gives a larger estimate than 
MDL, and modelers usually choose an order no smaller than 
the AIC prediction. Therefore, an order of 15 has been chosen 
in [40] to extract the wave features of interest for both metallic 
and coated cones. 











B. Estimation of Signal-to-Noise Ratio 


The SNR of the radar system is influenced by noise 
contributed by several sub-systems, such as oscillator phase 
noise, mixer //f noise, thermal noise, antenna leakage and 
mixer leakage. It is cumbersome and often inaccurate to 
estimate these random noise levels in the radar sub-systems 
using circuit models [70]. SSM considers this noise through 
decomposition of the Hankel matrix into signal and noise 
subspaces using the singular value decomposition. 


The threshold for such decomposition is based on the 
dynamic range and SNR, which can be estimated from power 
spectrum of the original measured or simulated data, as 
discussed next. First, we consider measured data on the chest 
wall displacement of a human subject to extract the heart and 
respiration rates [41]. UWB radar with a center frequency of 
2.4 GHz, bandwidth of 2 GHz and pulse repetition frequency 
(PRF) of 75 Hz is used to collect the data. The complex signal 
I+jQ of the radar return from a human subject located 1 m 


from the radar is plotted in Fig. 2 [41]. Random fluctuations 
caused by measurement noise are evident. The complex 
signal in Fig. 2 1s transformed with FFT using a Hamming 
window to suppress the sidelobes. This does not affect the 


peak signal but clearly defines the noise floor relative to the 
peak. Therefore, the dynamic range and SNR can be estimated 
from the compressed (or transformed) signal plotted as a 
function of Doppler frequency proportional to the chest wall 
displacement. Fig. 3 shows that the SNR for Channel 3 data is 
moderately high (around 24 dB), as one may expect from 
indoor laboratory RCS measurements, which are not subject 
to multipath and clutter encountered by fielded radar systems. 
Next, we consider SNR evaluation for simulated data. 
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Fig. 2. Raw time sequences of I and Q channels [41 ]. 
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Fig. 3. SNR estimate from UWB radar data on a human subject [41]. 


In simulated data, the difference between the desired and 
the modeled signal, as defined in (1), is a combination of 
model mismatch error and any noise intentionally added to the 
signal to determine statistical effects of random variation in 
parameters (e.g., Monte Carlo trials). In the absence of any 
added noise, the model mismatch error can be treated as noise. 
The SNR can be defined as the ratio of variance of the signal 
data sequence to the variance of the noise sequence, 


oe } k=1,2,--,N (42) 
V i y(k)— W(k)} 


where V stands for variance defined by 


Vf} =——Y | f-H 
N-1- 


SNR = 1010| 


2 
’ 





(43) 


and mw is the mean of f(k), given by 
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1 N 
w= YS) (44) 


In a noise-free situation, if the SSM estimates defined by 
(34) closely match the data samples, the SNR defines the 
dynamic range of the model. Thus, the number of scatterers or 
signals giving rise to the estimates corresponds exactly to the 
number of sinusoids embedded in the data. In a noisy 
situation, it is important to use caution so that the estimates 
are not direct replicas of the noisy samples. Therefore, an 
“optimal” model order should provide appropriate SNR or 
dynamic range for a sinusoid being modelled. 


IV. RESULTS AND DISCUSSION 


In this section, two examples based on simulated data are 
presented to illustrate the state space method. In the first 
example, we consider a planar dielectric slab illuminated by a 
plane wave at normal incidence and isolate the reflection off 
the front face using range-classified poles pertinent to the 
specular reflection. The subsequent radar-return signals 
enable classification of multiple internal reflections also. In 
the second example, we consider Mie scattering from a sphere 
in the absence of noise, extract the creeping wave that 
circumnavigates the shadow zone and returns to the radar, and 
examine accuracy of the SSM estimates by comparison with 
an analytically derived expression of the creeping wave. Next, 
random white Gaussian noise is added to the Mie series 
solution, and Monte Carlo simulation is performed to examine 
robustness of the estimates to noise. Numerical 
considerations such as SNR, singular value spectrum and 
order determination are addressed in detail for both of these 
examples. The reader is referred to [40]-[43] for application 
of SSM to feature identification using monostatic RCS 
measured data on stationary targets as well as human subjects. 


A. Reflection by a Dielectric Slab 


Plane wave reflection from a 15-mm thick lossy dielectric 
slab with dielectric constant ¢, = (5,—0.01) is considered, and 
it is shown that the response of the front face of the slab as 
well as multiple internal reflections can all be isolated from 
the composite response using SSM. Each of these 
contributions is specifically mapped to a range-gated pole. 
The slab of thickness d is located in free space. Its Fresnel 
reflection coefficient is calculated as 


(n? -n tanh (7d) 
2MM +n +75 tanh (y,d) | 
7 =jk, k =kg8E,, ko = Oy 4h80., 


7), =F dE My =120z. 


Without loss of generality, normal incidence is considered for 
illustration. The complex propagation constant and intrinsic 
impedance in the lossy dielectric are given by y, and 7,, 


ya) =T= (45) 


(46) 


respectively. The phase constant and the intrinsic impedance 
of free space are denoted as k, and 7,, respectively, while œ 
is the angular frequency. The reflection coefficient is 
calculated at normal incidence over 2-20 GHz bandwidth and 
Fourier transformed, using a Hamming window to suppress 


the sidelobes, to obtain the composite response shown in Fig. 
4. The transform employs distance to the scatterer (or range) 


instead of time, an operation known as pulse compression. 
i 
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Fig. 4. Range response of the 15-mm thick slab reflection coefficient. 


Because the phase reference is at the front face, the pulse 
with peak at zero range 1s the reflection off the slab-front, the 
second pulse is the first reflection off the back face arriving 
coherently at the “receiver,” the third pulse is the second re- 
reflection off the back, and so on, as illustrated geometrically 
in Fig. 5 using transmission line analogy [71]. 
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Fig. 5. Reflection diagram in analogy with transmission line circuit 
[71]. 


Next, a 10" order state-space model (see (1)) is estimated 
from the simulated data. The model parameters are listed in 
Table 1, where the poles and the signal amplitudes are 
classified in terms of scatterers of interest highlighted in the 
last column. The rows are sorted in descending order of the 
peak amplitude. All the poles lie within or on the unit circle, 
validating the stability of the system. 


Table 1. SSM parameters for model order M = 10. 


Index Pole Location Complex Amplitude Abs(Amp) Range (m) Notes 


1.0000 - j0.0000 
1.0000 -j0.0070 
0.9999 -j0.0140 
0.9998 -j0.0211 
0.9996 -j0.0281 


0.9993 -j0.0351 
0.9994 -j0.0421 
0.9987 -j0.0492 
0.9984 -j0.0562 
0.9979 -j0.0632 


-0.3820 + j0.0004 
-0.3077 -j0.1057 


0.0374 + j0.0290 


-0.0038 - j0.0058 


0.0002 +j0.0010 
0.0000 - j0.0001 
0.0000 + j0.0000 
0.0000 - j0.0000 
0.0000 +j0.0000 
0.0000 -j0.0000 


0.381967 
0.325323 
0.047331 
0.006886 
0.001002 
0.000146 
0.000021 
0.000003 
0.000000 
0.000000 


0.0000 
0.0150 
0.0300 
0.0450 
0.0600 
0.0750 
0.0900 
0.1050 
0.1200 
0.1350 





Front 
Refl. 1 
Refl. 2 
Refl. 3 
Refl. 4 
Refl. 5 


It is interesting that the signal peaks are separated in SSM- 
computed range (see (38)) by exactly 15 mm, the thickness of 
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the slab. The dominant pole, located at (1, 0) with amplitude 
of 0.382 and phase z, is found to have zero range. This signal 


corresponds to the isolated front face reflection. The 
subsequent poles, consecutively spaced at 15 mm, represent 
the isolated secondary reflections off the back face, as 
illustrated in Fig. 4. The first five signals of interest in Table 
I appear to have the maximum spectral content, as the sixth 
term is three orders of magnitude smaller than the dominant 
signal. Using the coherent summation of only these first five 
terms in the state space model (1), we have computed the 
frequency response of the composite signal. Fig. 6 compares 
the “truth” with the state space model in both magnitude and 
phase, and the two sets of curves overlap each other, 
signifying excellent model fidelity. 


As alluded, one can isolate the front face response by using 
only the dominant pole at zero range. The reflection 
coefficient of this isolated pulse is the same as that of the half- 
space problem because the delayed reflections off the back 
face are not included. Therefore, using the constant amplitude 
of I’, =—0.382 + 70.004 from Table 1, the slab’s dielectric 


constant can be determined as 
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(b) 
Fig. 6. Corroboration between SSM model and the truth (eq. (45)) 
using only the first five signals in Table 1. (a) Magnitude, (b) Phase. 
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£, = L|; (47) 
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Eq. (47) yields £. = 5-— j0.01 exactly, verifying the accuracy 


of the signal extraction. It is imperative that isolation of 
individual signals of interest requires adequate bandwidth to 
provide the spectral resolution needed to separate the peaks. 
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Fig. 7. Range response of the 5-mm thick slab reflection coefficient. 
When the thickness of the slab is reduced to 5 mm, one 
obtains the compressed range response shown in Fig. 7, which 
clearly shows the impending merging of the first two peaks. 
Using a bandwidth of 18 GHz, these two peaks are very close 
in amplitude, and may not be adequately resolved in range. 
But when the bandwidth is increased to 28 GHz, the peaks can 
be isolated. Furthermore, in the latter case, the markers also 
indicate secondary fluctuations appearing at approximately 10 
and 15 mm, corresponding to re-reflections off the back face. 
However, in measured data, such small fluctuations of the 
main lobe may be masked by noise, and a larger bandwidth 

may not necessarily improve the SNR. 


Fig. 8 plots the model error magnitude of the frequency- 
domain response for SSM with order M = 5 (the first five 
terms in Table 1) and d= 15 mm. The mean error is -77.3 dB, 
proving that all the significant signals have been modeled. A 
similar agreement between raw data and SSM is observed in 
the phase too, with mean error of 0.2°. 


Frequency Domain Error, Abs(Raw Data - SSM), d= 15 mm, M=5 
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Fig. 8. Model error in magnitude of the SSM frequency response. 


The model error for M = 10 and M = 15 did not decrease 
beyond the error depicted in Fig. 8 for M = 5. A performance 
measure for evaluation of the model error is the SNR in (42). 
Table 2 lists the SNR for few model orders ranging from 5 to 
15. Itis evident that SSM with M = 10 gives the highest SNR 
and any further increase in the order does not impact the error. 


Table 2. SNR as a function of model order for the slab problem. 


A deeper understanding of the sources of modeling error in 
SSM can be obtained by examining the error for various 
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model orders relative to the compressed SSM signal. Along 
with the singular value spectrum, to be discussed shortly, this 
will help us in evaluating impact of the signals discarded in 
the model as “noise.” Fig. 9 displays the error responses for 
model orders M = 5 and M = 10, compressed using a 4,096- 
point FFT on sequences of length 3,601, with a Hamming 
window to suppress the noise sidelobes. For comparison, the 
compressed SSM signal estimate, identical for the two model 
orders, is also shown. The range coordinate for each peak (in 
mm) is annotated in the graph. We observe that the SSM 
signal, (k), models the five dominant reflections signified by 


the first five poles (see Table 1), and the remaining signals are 
embedded in noise. Because these discarded signals are there 
in the truth data (45), it is not surprising that the range spectra 
of the error signal, y(k)— (kK), depict these discarded noise 


peaks precisely. The number of such noise peaks embedded 
in the error sequence depends on the model order. For 
example, with M = 5, we calculate a noise floor of -110 dB 
(the mean square of the compressed error pulse) and identify 
the three noise peaks at 75, 90 and 105 mm. With M = 10, we 
can identify the noise peaks at 105, 120 and 135 mm. 
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Fig. 9. Comparison of the compressed signal and noise (error) 
responses for two model orders. The numbers above each peak 
indicate the range in mm. 


The noise spectrum did not change with any further increase 
in model order beyond 10. Thus, the last noise peak one can 
identify is at 135 mm. As the highest noise peak (at 75 mm) is 
90 dB down from the main reflection at zero range, none of 
these discarded noise peaks has any influence on the SNR. 
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Fig. 10. Singular value spectrum to determine the model order. 


Lastly, we estimate the model order by examining the 
singular value spectrum plotted in Fig. 10 as a function of its 


index, r. The SVs are normalized to the largest value (r = 1). 


The model order can be determined subject to the threshold 
[56] o,/o,+10”, where p is the number of significant 


decimal places in the truth data. Thus, if the data is accurate 
to 7 places, the SV spectrum predicts a model order of 7 = 10. 
In practice, the data may be accurate to third or fourth decimal 
place, but not till the seventh as the simulated reflection 
coefficient (45) for this simple illustrative problem. For three- 
digit accuracy, Table 1 indicates that a model order of 10 
yields 5 signal poles and 5 noise poles (weak signals 
embedded in noise). Indeed, the error analysis in Figs. 8 and 
9 validates this observation. 


B. Scattering by a Sphere 


Let us consider the extraction of creeping waves from the 
canonical problem of scattering by a perfect electrically 
conducting (PEC) sphere using the SSM. The scattered field 
for plane wave incidence is computed analytically using the 
Mie series [72], first in the absence of noise, and then with the 
influence of additive Gaussian noise. The model order is 
estimated in each case by examining the singular value 
spectrum, and the accuracy of the extracted creeping wave is 
established by comparison with an analytically derived 
asymptotic expression of the creeping wave [73]. 


B.1. Baseline without Noise 


The monostatic RCS is calculated for a PEC sphere of 
radius a =17.7 cm using the Mie series, and is plotted against 
frequency over arange of 2 to 20 GHz in Fig. 11. As the radius 
becomes much larger than the wavelength, the normalized 
RCS asymptotically approaches ma. The contribution to 
back-scattered field comes from a specularly reflected ray 
emanating from point A, and a creeping wave, which attaches 
to the sphere tangentially at point B, navigates half the 
circumference, and detaches at the symmetrical tangential 
location C, as shown in the inset in Fig. 11. The specular 
contribution is governed by geometrical optics (GQ) 
approximation and can be analytically calculated [74]. The 
creeping wave diffracted field can be approximated by 
asymptotic evaluation of the Fock integrals [73], [74]. 
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Fig. 11. Normalized RCS of a PEC sphere, computed using Mie 
series. 
In order to isolate the scattering centers pertinent to the 
specular and the creeping wave, we first compute the range 
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profile of the back-scattered field, displayed in Fig. 12, by 
FFT of the Mie series frequency response using a Hamming 
window to suppress the sidelobes. The phase reference is 
chosen such that the point of specular reflection is at zero 
range. The second peak which corresponds to the creeping 
wave is at 0.454 m range, in agreement with the range 
predicted from ray path geometry (see the inset of Fig. 11) as 
0.455 m, i.e., R =(7z/2+1)a. 
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Fig. 12. Range profile of the back-scattered field. 


Next, SSM is applied to extract range-isolated poles 
specific to the creeping wave, and their cumulative frequency 
response is calculated. A 10" order state space model is 
computed from the Mie series frequency response shown in 
Fig. 11. We have observed excellent corroboration between 
the estimated SSM model and the truth data in Fig. 11 over 
the entire 18 GHz band. For brevity, we focus only on 
representation of the creeping wave using the corresponding 
range-isolated poles from the SSM. As described in Section 
II, an advantage of using range processing with the spectral 
estimation method is the direct relation between range and 
pole phase (see (38)). Thus, one may isolate a given scattering 
mode by adding contributions from only the poles associated 
with the range window of that mode. Of the 10 poles used in 
the SSM to represent the entire signal in Fig. 11, only two 
poles are identified with model-computed range of 0.454 m, 
relevant to the creeping wave peak. By coherently summing 
the contributions of only these two poles, we obtain the 
extracted frequency response shown in Fig. 13. 
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Fig. 13. Comparison of SSM-extracted creeping wave with the 


analytical solution from [73]. 


In order to demonstrate baseline validation with no noise 
present in the data, we also plot in Fig. 13 an analytical 
solution for the spherical creeping wave from [73], which 
essentially overlaps with the SSM-estimated data. The signal 


amplitudes for the two creeping wave poles are observed to be 
around -36 dB, signifying high accuracy in the model even for 
small signals. Later, we will evaluate model robustness using 
Monte Carlo analysis with Gaussian noise added to the Mie 
series. The absolute model prediction error for the creeping 
wave estimation is plotted in Fig. 14. A similar agreement 
with the analytical solution is also observed for the specular 
peak extracted from the Mie series in Fig. 11. 
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Fig. 14. Model prediction error for the extracted creeping wave 


relative to the analytical solution from [73]. 
B.2. Monte Carlo Analysis 


The SSM extraction of signals of interest from Mie series 
in the presence of noise is considered next. Monte Carlo (MC) 
analysis is performed by adding random noise of a given SNR 
to the Mie series solution for the PEC sphere and evaluating 
the accuracy of the SSM-extracted components, namely, the 
creeping wave and the specular. The measurement noise w(k) 
(see (1)) is assumed to be complex white Gaussian with 


variance o7, defined by peak signal-to-noise ratio 


(48) 


O 


n 


SNR = 20 e| Z= | 


where o}, denotes the signal variance. For a given SNR, 1000 


independent trials are executed on the Mie series solution, and 
the state space method is applied to process each noise- 
corrupted data sequence and extract the wave constituents of 
interest. In order to assess the quality of the data, the noise- 
corrupted Mie series (truth) data is plotted in Fig. 15 against 
SSM model (with order M = 10) of the total response for the 
mean of 1000 MC trials with SNR = 20 dB. Itis observed that 
coherent averaging of the MC trials reduces the influence of 
noise considerably. As the electrical size of the sphere 
increases, the asymptotic limit of SSM correctly reaches za’. 
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Fig. 15. Average of the SSM-extracted Mie series composite signal 

over 1000 Monte Carlo trials. 
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Fig. 16. Average of the SSM-extracted specular wave signal over 
1000 Monte Carlo trials. 


Next, we consider the extraction of the specular wave. Fig. 
16 depicts the frequency response corresponding to SSM 
processing of the average of 1000 MC trials, with SNR 
ranging from 5 dB to 20 dB. Only two poles are used in the 
SSM for the extraction. The response for each SNR is 
observed to track the analytical (noise-free) solution quite 
well. The worst-case error relative to the reference solution is 
about 0.05 dB and occurs for SNR = 5 dB. More importantly, 
because of the relatively large amplitude of the specular 
(about -10 dB), we have observed good correlation between 
corresponding pole locations for each trial. For brevity, the 
pole plots are not included. 


We address Monte Carlo analysis of the creeping wave 
next. As seen in Fig. 13, the noise-free creeping wave signal 
has an amplitude of -20 to -70 dB, and therefore, adding noise 
before the creeping wave extraction would considerably stress 
the state space algorithm. Fig. 17 displays the performance of 
SSM in extracting the creeping wave for various SNRs, 
relative to the reference analytical solution from [73]. 
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Fig. 17. Comparison of SSM-extracted creeping wave for the 
average of 1000 MC trials with the analytical solution from [73]. 

The extracted signal is in reasonable agreement with the 
analytical solution for SNR 210 dB. The SSM estimate 
deviates significantly from the reference solution at 
frequencies greater than 10 GHz for SVR = 5 dB. For SNRs 
between 5 and 15 dB, the error increases with frequency for 
f >16 GHz, but the RCS response is around -70 to -65 dB, 


which is approaching the noise floor. It has been observed that 
the pole drift from trial to trial becomes significant for SVR = 
5 dB, and in fact, for signals smaller than -60 dB, it becomes 


difficult to discriminate noise poles from the signal poles. 
Therefore, for weak signals such as creeping waves, caution 
should be exercised in extracting the signal under low SNR 
conditions. Nevertheless, the worst-case performance of SSM 
under the stressing conditions depicted in Fig. 17 is gratifying, 
given that the frequency response of the creeping wave signal 
for each SNR, extracted from the average of 1000 MC trials 
of the Mie series, follows the overall analytical trend. One can 
improve the noise performance by employing two- 
dimensional data, e.g., aspect- and frequency-dependent RCS, 
which improves the SNR by coherent integration of frequency 
samples over many pulses [61], [62]. Our extensive work on 
state-space methods to process biomedical radar data [41 ]- 
[43], [75], [76] suggests that the SNR of vital sign detection 
can be substantially improved, and additional features such as 
subject localization, motion compensation, gait analysis, etc., 
can be estimated, using block-processing of one-dimensional 
data (in the frequency domain) to attain considerably 
improved accuracy over FFT-based methods. 
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Fig. 18. Model prediction error for the extracted creeping wave 
relative to the analytical solution from [73] after 1000 MC trials. 

The absolute model prediction error for the creeping wave 
estimation from the average of 1000 MC trials on the Mie 
series is plotted in Fig. 18 in terms of the SNR. It is evident 
that the worst-case error for SVR = 5 dB is between two to 
four orders of magnitude larger than the model error for the 
noise-free case shown in Fig. 14. For SVR =10 dB, the error 
is only slightly larger than that in Fig. 14. 
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Fig. 19. Singular value spectrum vs. model order for noisy data with 
SNR = 5. MDL yields a model order of 7 and AIC, an order of 20. 


Lastly, we investigate the model order by computing the 
singular value spectrum as well as AIC and MDL estimates. 
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It is seen in Fig. 19 that the transition from strong signal 
components to small singular values is not well-defined 
because of noise. AIC yields model order of 20 while MDL 
yields 7. Due to measurement noise, AIC usually gives a 
larger estimate than MDL, and modelers usually choose an 
order no smaller than the AIC prediction. Therefore, an order 
of 20 has been chosen to extract the wave features of interest 
for the sphere. Compared to estimates in the dielectric slab 
case, which have a very high SNR (see Table 2), the accuracy 
in the SSM estimate for the MC trials is much lower, of the 
order 0.001. 


V. CONCLUSIONS 


A spectral model based on state-space ARMA 
representation has been presented to isolate and extract modal 
EM responses, such as creeping waves and multiply reflected 
or diffracted waves, which are of interest in radar target 
identification and feature extraction. The eigenvalues of the 
open-loop system matrix, i.e., poles of the ARMA transfer 
function, provide the range and the frequency-dependent 
amplitude decay/growth rate of the field data. The amplitudes 
in the complex exponential model are obtained by least 
squares modal decomposition involving the state matrices. 
Range classification of the poles allows for isolation of 
desired responses in the EM signature via pulse compression 
and spectral decomposition. For purposes of illustration, in 
this paper the SSM is reviewed by application to two simple 
EM problems, namely, Fresnel reflection by a dielectric slab 
and the extraction of creeping waves using Mie scattering by 
a PEC sphere. The method has been applied to isolate the 
leading-edge reflection and multiple internal reflections of the 
slab. It is shown that range-classification of the singularities 
in the response enables characterization of the isolated 
reflected waves using a low model-order SSM representation. 
In the second example, we have analyzed the extraction of 
creeping waves from Mie scattering by a PEC sphere in the 
absence of noise. The SSM estimates of the extracted 
creeping wave have been validated by comparison with an 
analytically derived expression of the creeping wave. Next, 
random white Gaussian noise is added to the Mie series 
solution, and Monte Carlo simulation is performed to examine 
robustness of the SSM estimates to noise. Numerical 
considerations such as SNR, singular value spectrum and 
order determination are addressed in detail for both of these 
examples. The reader is referred to [40]-[43] for application 
of SSM to feature identification using monostatic RCS 
measured data on stationary targets as well as human subjects. 
As we have shown in [40], the isolation of electromagnetic 
wave species of interest, such as creeping waves, multiply 
diffracted waves and specular scattering, yields a better 
understanding of the physics behind wave propagation around 
curved dielectric or coated structures, thereby improving the 
accuracy of feature extraction or target identification. 


APPENDIX 


In terms of the  finite-rank observability§ matrix 
Q e COD computed from the SVD in (25) as 


Gale ca CP = CA cH]. (49) 


the state transition or the open-loop matrix AeC"”™ is 
determined by the solution to the matrix equation 


GAO... (50) 
with 
Õ „=[C4 Ce CH rey. 6D 
= 2 were! 
O.,,=|C CA CA o | (52) 


It is observed that the matrices Q_, and Q_,, are obtained by 


deleting the first and last rows, respectively, of the matrix Q 
in (49). Using the pseudo-inverse least squares on (50), we 
obtain the state transition matrix A given in (26) and repeated 
below. 
“e x -1 x, x 
A= (aa) GF, 2... (53) 
This matrix may also be derived from the controllability 
matrix [I c C% derived from the SVD in (25) as 
[=|B AB AB AB ANB, (54) 


It follows that the state transition matrix A satisfies the matrix 
equation 


AL eT (55) 

with 
a =|4B ÆB 4B AB], (56) 
Ta =|B AB AB A™°B |. (57) 


Eqs. (56) and (57) are obtained by deleting the first and last 
columns, respectively, of the controllability matrix I in (54). 
By solving (55) for A using the pseudo-inverse least squares 
method, equivalent to (53), we obtain (also see (27)) 


A= ks (agg r i 


(58) 


Next, we address the computation of control and 
observation matrices, B and C, respectively, using two 
alternative approaches. In the first approach, it follows from 
(49) that the observation matrix C is simply given by the first 
row of the observability matrix, 


C= Qaia; (59) 


Alternatively, the IIR approximation, (k), of the data 
sequence y(k), k=1,2,---N, can be computed in terms of 
the state matrices using (16) as 

p(k) =CA*'B, k=1,2,---N. (60) 


We define the augmented observability matrix employing all 
the NV samples: 


Õ,=[C CA CH uE (61) 
which is related to the impulse response in (60) as 
Q B=F. (62) 
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Therefore, in order to minimize the error between the model 
y and the measured data vector y, the control matrix B is 


computed by the pseudo-inverse least squares method, as 
shown below. 


B=(O10,) O43". (63) 


In summary, in the first approach, one may calculate the 
matrices C and B from the observability matrix Q, using (59) 
and (63), respectively. 

Alternatively, in the second approach, B is computed from 


the first column of the controllability matrix in (54) as 


B=T(:,)). (64) 


and C follows from the least squares fit between the data 
sequence y(k), and the state-space IIR approximation y(k) 


in (60), as shown next. In terms of the augmented 
controllability matrix 


D, = |B AB AB ANB|, (65) 
the estimation problem for C may be written as 
CT y= 3, (66) 
and its pseudo-inverse least squares solution yields 
CSG... (67) 


To summarize, it 1s emphasized that there are two 
expressions for the state transition matrix, A, and two 
corresponding alternative approaches to computing B and C. 
The matrix A may be computed from either the observability 
matrix using (53), or the controllability matrix using (58). 
Correspondingly, the matrices B and C must be computed by 
using either (63) and (59), or (64) and (67), respectively. In 
order to improve the estimation accuracy in a low SNR 
environment, we almost always employ the least-squares 
computations (63) and (67) to estimate B and C, respectively, 
and seldom use (64) and (59). 
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